#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>

// Número de iterações = 10^9
#define N 1000000000
// Número de threads
#define N_THREADS 8

pthread_t threads[N_THREADS];
pthread_mutex_t mutex_sum;
int p = 0;
unsigned short xsubi[3];

// Função que cada thread acessa para calcular o PI
void *calculatePI(){
	double x, y;
	int i, local_p = 0;
	int i_final = N/N_THREADS;

	// Como são utilizadas N_THREADS, o intervalo do loop pode ser então dividido entre elas N/N_THREADS
	for(i = 0; i < i_final; i++){

		// Geração aleatória dos pontos (x,y)
		x = erand48(xsubi);
		y = erand48(xsubi);

		// Equação da circunferência de raio 1 (verifica se [x,y] estão dentro dela)
		if((x*x) + (y*y) <= 1.f){

			// Para que não seja necessário utilizar o lock do mutex "mutex_sum" em todas as iterações do loop para cada thread, é feita uma soma parcial
			// Dessa forma, como as variáveis locais não são compartilhadas entre threads, só é necessário um lock para atualizar a soma final
			local_p++;
		}
	}
	
	// Para que somente uma thread accesse a variável global "p", é utilizado o mutex "mutex_sum". Desta forma, a variável "p" é atualizada corretamente
	pthread_mutex_lock(&mutex_sum);
	p += local_p;
	pthread_mutex_unlock(&mutex_sum);

	pthread_exit(NULL);	
}

int main(int argc, char *argv[]){

	double pi;
	int t, rc;

	pthread_mutex_init(&mutex_sum, NULL);

	for(t = 0; t < N_THREADS; t++){
		// Criação das threads
		rc = pthread_create(&threads[t], NULL, calculatePI, NULL);
			
		if(rc){
			printf("ERROR: %d\n", rc);
			return -1;
		}
	}

	// Junção das threads quando completarem sua execução (em ordem de criação)
	for(t = 0; t < N_THREADS; t++)
		pthread_join(threads[t], NULL);

	// Cálculo do PI
	pi = 4.f * (((double) p) / ((double) N));
	printf("%.6f", pi);

	pthread_mutex_destroy(&mutex_sum);
	return 0;
}
